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We propose an extension of the numerical approach for integrable Richardson-Gaudin models 
based on a new set of eigenvalue-based variables 1 2 . Starting solely from the Gaudin algebra, the 
approach is generalized towards the full class of XXZ Richardson-Gaudin models. This allows for 
a fast and robust numerical determination of the spectral properties of these models, avoiding the 
singularities usually arising at the so-called singular points. We also provide different determinant 
expressions for the normalization of the Bethe Ansatz states and form factors of local spin operators, 
opening up possibilities for the study of larger systems, both integrable and non-integrable. These 
expressions can be written in terms of the new set of variables and generalize the results previously 
obtained for rational Richardson-Gaudin modelJ^ and Dic.ke-Jaynes-Cummings-Gaudin models 4 . 
Remarkably, these results are independent of the explicit parametrization of the Gaudin algebra, 
exposing a universality in the properties of Richardson-Gaudin integrable systems deeply linked to 
the underlying algebraic structure. 

PACS numbers: 02.30.Ik,02.60.Cb,21.60.Fw,74.20.Rp 


I. INTRODUCTION 

Exactly-solvable systems play an important role for 
the understanding of the physics of quantum many-body 
systems. They offer an insight into the behaviour of 
strongly correlated systems in ways that would otherwise 
be impossible. One class of these systems is the class of 
Richardson-Gaudin (RG) integrable systems, which can 
be derived from a generalized Gaudin algebra 5 The 
pairing model in the reduced Bardeen-Cooper-Schrieffer 
(BCS) approximation, used to describe superconductiv¬ 
ity, has been shown to be RG integrable^, as has the p x + 
ip y pairing HamiltoniarP®, the central spin modeP®. fac- 
torisable pairing models in heavy nuclei tH, an extended 
d + id pairing Hamiltoniarff 2 and several atom-molecule 
Hamiltonians such as the inhomogeneous Dicke modeF®. 
For these models, diagonalizing the Hamiltonian in an 
exponentially scaling Hilbert space can be reduced to 
solving a set of non-linear equations scaling linearly with 
system size. For most of these systems, computationally 
inexpensive expressions are also known for several form 
factors and overlaps, which can be used to investigate ob¬ 
servables in systems where the Hilbert space is too large 
for traditional exact methods. Unfortu nately, the Bethe 
Ansatz or Richardson-Gaudin equations 10 1 iffi]that need 
to be solved are highly non-linear and give rise to sin¬ 
gularities, maki ng a straightforward numerical solution 
challenging 16 17 [ Several methods have been introduced 
as a wa y to resolve this difficulty, such as a chang e of 
variables 17 EH, a (pseudo-)deformation of the algebraP^Sl, 
or a Heine-Stieltjes connection, reducing the problem to 
a differential equation 21 . The ground state energy in the 
thermodynamic limit has also been obtained by treating 
the interaction as an effective temperature!^®. 

Our method consists of a generalization of the numeri¬ 
cal method first proposed by Faribault et all® for a set of 


non-degenerate rational RG models, and later extended 
to degenerate rational models by El Araby, Gritsev and 
FaribaulfP. The equations for the Dicke model were in¬ 
dependently presented by Babelon and Talalaev 23 ! In 
this method, an alternative set of equations is derived in 
terms of variables related to the eigenvalues of the con¬ 
stants of motion. One of the advantages of this method is 
that it offers not only the solutions of the RG equations, 
but also effi cien t numerical expressions for form factors 
and overlaps 3 ®. Such expressions for the rational model 
have already been used to study quantum quenche^Sl and 
decoherence in quantum dotiP® anc j should facilitate the 
use of the XXZ model as a variational or projected ap¬ 
proach for non-integrable systems, similar to the use of 
the XXX model in quantum chemistrji^®^® 

In the following, the results for the rational XXX model 
are generalized towards the full class of XXZ RG inte¬ 
grable systems and possible applications are discussed. 


II. RICHARDSON-GAUDIN MODELS 

A. Definitions 

Unlike its classical counterpart, no single unique def¬ 
inition is known for quantum integrability 29 . One class 
of systems usually denoted integrable is the set of RG 
integrable systems, since these support as many (non¬ 
trivial) conserved operators commuting with the Hamil¬ 
tonian as there are degrees of freedom in the systeuPS®. 
These systems can be diagonalized exactly by means of a 
Bethe Ansatz wave function, where diagonalizing the full 
Hamiltonian can be reduced to solving a set of nonlinear 
coupled equations for the variables in the ansat^EE®. 

The families of RG integrable syst ems have their roots 
in a generalized Gaudin algebra® 1 based on the su( 2) 
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algebra of (quasi-)spin operator^. Depending on the 
physical systems, they can either represent spin degrees 
of freedom, fermion pairs or bosonic/fermionic Schwinger 
representations. The relevant (quasi-)spin operators can 
be defined as satisfying the following commutation rela¬ 
tions 

[S?,lS\] = 6 ij sl [S?,S j ] = -6 ij S i , 

[S},S J } = 2S. lJ S^, (1) 

where each separate algebra spans a su(2) i algebra asso¬ 
ciated with a level i and irreducible representations (ir- 
rep) | di,Hi). For a set of n levels, the RG integrable 
models are then defined by a set of n mutually commut¬ 
ing operators 


Ri — S i + g 




-x ik (sjs k + SiSD + z ik s"s u k 


( 2 ) 


Following Gaudiifl® and Dukelsky et alP-^, it is possible 
to obtain a set of conditions for which the set of operators 
Ri commute mutually. These are given by 

Xij — Xj i . Zij Zjj . (3) 

XijXjk — Xi k {Zij + Zj k ) = 0, (4) 

the so-called Gaudin equations. They were originally dis¬ 
covered by Gaudin as defining a general class of quadratic 
Hamiltonians in the spin variables, among which the 
Gaudin magnet. The derivation by Dukelsky et al. dif¬ 
fers from the derivation by Gaudin in the presence of the 
generator of the Cartan subalgebra S' 0 , which however 
does not influence these conditions. In fact, the Gaudin 
models are given by Dukelsky’s constants of motion in 
the limit g —> oo . 

Gaudin mentioned three classes of solutions of these 
equations, where all classes consider X{j and Zij as 
odd functions of some real arbitrary parameters Xij = 
X(ei,€j). The physical interpretation of these parame¬ 
ters follows from the expression for the Hamiltonian con¬ 
structed with these parameters. 


1. The rational model 


Xij Zij 



( 5 ) 


2. The trigonometric model 

v _ 1 

13 sin (ej-ej)’ 


Z^ = cot(e,; - ej) (6) 


3. The hyperbolic model 

1 

13 sinh(ej — ej)' 


Z^ = coth(e; - ej) (7) 


Here the rational model is also referred to as the XXX- 
modef^, indicating that the coefficients in the expression 
for the constants of motion are identical for all 3 com¬ 
ponents of the spin. In the same vein, the trigonometric 


and hyperbolic models are called XXZ-models. Alterna¬ 
tive solutions are given by 


Xij — 


Zij — 


\/l + 2 u€i + fief 1 + 2 aej + fief 


1 + a(ei + £j ) + , 


c * c j 


( 8 ) 


which was proposed by Richardson^, and 


x„ = Z„ = (9) 


which is a reparametrization of the hyperbolic modeP. 


B. Diagonalizing integrable Hamiltonians 


Any general RG integrable Hamiltonian can be written 
as a linear combination of the constants of motion 

n 

H = Y,ViRi, (io) 

i=1 


resulting in a wide variety of systems. The reduced BCS 
Hamiltonian can be found from the rational modeP, the 
Px + Wy pairing Ha miltonian can be derived from the 
hyperbolic modeP^, and the central spin model Hamil¬ 
tonian is identical to one of the constants of motion of 
the rational model 10 . By introducing a bosonic degree of 
freedom, the inhomogeneous Dicke model can be found 
as a limiting case of the XXZ modefP 

Since all constants of motion commute mutually, they 
also commute with any Hamiltonian that can be writ¬ 
ten as Eq. (10), so this reduces the cliagonalization of 
the Hamiltonian to the diagonalization of one of the con¬ 
stants of motion. It can be shown that the eigenstates 
are given by a Bethe Ansatz 


n / n 

=n 

a=l \i =1 

with N the number of excitations^ and the vacuum state 
1 9) = g>" =1 \di, — di) the tensor product of the lowest- 
weight irreps of each su(2)i copy, provided the RG equa¬ 
tions 



N 


1 T g ^ ) Zi a di g ^ ) Zg a — 0, 

2—1 /3^a 


Va = 1 ... TV, (12) 


are satisfied, where the algebra has been extended to 
Xi a = X(ei,x a ) and Z ia = Z(ei,x a ) by introducing 
a set of TV parameters {x a }. These parameters, also 
known as the RG variables or rapidities, fix the wave 
function and need to be determined from the RG equa¬ 
tions (12). There are multiple strategies to derive this 


result. Richardson obtained these equations for the re¬ 
duced level-independent BCS model starting from a vari¬ 
ational approact l 14 F 4 ( Gaudin started from integrability 
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constraints-® Zhou et al. derived these results starting 
from the Algebraic Bethe Ansat^ 5 and Ortiz et al. used 
a commutator scheme based on the Gaudin algebrgP. 

Once the RG equations have been solved and the ra¬ 
pidities determined, the eigenvalues of the constants of 
motion Ri can be evaluated as 


( n N \ 

1 T g 'y ) z ik d k + g y ) Zpi I , Vz — 1... . 
kjti / 3=1 / 

(13) 

Although the set of RG equations seems to be linear in 
the elements Zi a and Z a p, these variables are still cou¬ 
pled through the Gaudin equations, leading to a set of 
coupled and nonlinear equations. As an example, the RG 
equations for the doubly-degenerate (di = 1/2, Vi) XXX 
model are given by 



g Y ---=0, Va = 1... N. 

fr' X P ~ x <* 


(14) 

These equations need to be solved for the set of rapidities 
{x a } and are coupled and highly nonlinear. It has been 
shown that singular points arise in these equations at 
certain values of the coupling constant, where multiple 
rapidities x a equal one of the single-particle levels e LL£Hl£|. 
It can be readily seen that both the second term and the 
third term in the RG equations contain diverging poles, 
but it can also be shown that these singularities cancel 
exactly. Unfortunately, the behaviour arising at these so- 
called singular points hampers straightforward solutions 
of these equations, so involved numerical methods have 
to be found to circumvent these points. 


III. AN EIGENVALUE-BASED NUMERICAL 
METHOD 


A. Doubly degenerate models 


For the XXX spin-1/2 model (di = 1/2, Vz), it is pos¬ 
sible to introduce a new set of variables 

N 1 

Ai = A(ej) = ^2 -, (15) 

€i Xa 

a—1 


circumventing the singular points in the Richardson- 
Gaudin equation^ 1 . A set of equations equivalent to the 
RG equations can be found for these variables, however 
void of singular behaviour. As an illustration, the set of 
RG equations for the doubly-degenerate XXX model (Eq. 
14) are isomorphic to the set of quadratic equations 1 


A" = —A, 


V- Aj - Ai 




G - e, ; 


Vz = 1... n. 


(16) 


A straightforward numerical solution of these equations 
can easily be implemented and no singular behaviour will 
arise since no variables occur in the denominator. We 


wish to extend these results to the class of more general 
XXZ models, relying only on the Gaudin algebra and not 
on an explicit rational expression of A' and Z. Define 

N N 

A i = yz ia = yz(e i ,x a ), (17) 

a—1 a=l 

which reduces to the previously-introduced variables in 
the XXX model. It is interesting to note that the eigen¬ 
value ri of the constant of motion Ri is related to the 
variable Aj as 


n = di -1 - gAi + 


g Zikdk 

k^i 


(18) 


which has led to the denomination eigenvalue-based vari¬ 
ables. 

In the following, we will start from the Richardson- 
Gaudin equations for spin-1/2 systems 


n N 

1 + | y j Z ia = g y j Zp cn Va=l ...IV, (19) 

i—1 /3^CZ 

and several properties that can be derived from the 
Gaudin algebra, as previously obtained by Ortiz et alP. 
Firstly, it has been shown that 

(Xij) 2 — (Zij) 2 = T, Vz^j, (20) 

where T is a constant. The rational model corresponds 
to F = 0, while positive and negative F result in the 
trigonometric and hyperbolic models respectively. It also 
follows from the Gaudin algebra that 

ZijZjk — Z ik (Zij + Zjh) =r, (21) 

which is a generalization of the Gaudin equations for the 
rational model, where Xjj = Zjj and F = 0. Using only 
these equations, it can straightforwardly be shown that 


Aj = N(n — N)T - Ai + y Z,,{\j A,). V* = l...n 


9 


/ v -'J* 

3& 


( 22 ) 

with N the total number of excitations and n the number 
of single-particle levels. The full derivation can be found 
in Appendix [Aj Unlike the case of the rational model, 
these equations depend explicitly on the number of ex¬ 
citations N. When solving these equations numerically, 
it was found that the total number of solutions exceeds 
the dimension of the Hilbert space for N excitations dis¬ 
tributed over n levels. Therefore, these equations neces¬ 
sarily support unphysical solutions, not corresponding to 
any eigenstate, implying that this new set of equations 
is not equivalent to the original set of RG equations (Eq. 


121. In order to obtain a set of equations equivalent to 


the original equations, additional constraints for the Aj 
are needed. 

It is clear from the T = 0 case (XXX) that the new set 
of equations can not distinguish between the different 
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excitation sectors N. This can be imposed by noting 
that the sum of all constants of motion is given by the 
operator counting the number of excitations 

n n 

= ( 23 ) 

i=1 i=l 

Writing out the eigenvalues of the constants of motion in 
the new variables results in 

n 

- 9 -Y j K = n - (24) 


In the following section it will be shown that the full 
set of equations (22) and ([24]) has as many solutions 


as the dimension of the Hilbert space, so introducing 
this additional equation leads to a system of equations 
fully equivalent to the original set of RG equations. 
These eigenvalue-based equations do not show singular 
behaviour and can easily be solved numerically, as will 
be shown in the following. 


B. The weak-coupling limit 

In order to obtain some insight in the behaviour of 
the solutions of these equations, an approximate solution 
can be found in the weak-coupling limit (small |g|). This 
can be done by proposing a series expansion in g for the 
solutions of these equations and solving the equations at 
each order. For small g 7 a series expansion of A, in g can 
be proposed up to 0(g) 7 keeping only the two dominant 
terms 


d-d 


A,; = 


+ xf+0(g). 


(25) 


By plugging this expansion in the equations, we obtain 


+ 2 

+ O(g 0 ) = 0. 


+ i(2Af ) (l + A 


A^-Ar; 


This results in 


A- 1} = 0 or - 2 


and 


A® = 


2A.- 1] +2 


X>(W 

I#* 


- 1 ) _ a (-D 


(26) 

(27) 

(28) 


In order to satisfy Eq. ( |24| ), the number of dominant 
terms different from 0 has to equal the number of excita¬ 
tions N, resulting in (^) solutions. The total number of 
solutions then equals the dimension of the Hilbert space 
for N excitations distributed over n doubly degenerate 
levels. Any solution in the weak-coupling limit can be 


adiabatically connected to a solution for arbitrary cou¬ 
pling, indicating that all possible solutions are always 
found and no unphysical solutions are present. 

This series expansion can also be connected to the se¬ 
ries expansion obtained for the rapiditieJ®. In the limit 
g —» 0 we obtain a noninteracting model, where the ra¬ 
pidities x a converge to the parameters e^, depending on 
the corresponding distribution of excitations over energy 
levels. A rapidity converging to et then corresponds to an 
excited level i in the noninteracting limit. For finite but 
small g 7 dominant corrections of O(g) are present, which 
are proportional to the roots of orthogonal polynomials 
via a Heine-Stieltjes connectiorPSEZl. For x a converging 
to Cj, this results in Zi a = Z(ei 7 x a ) diverging as 1/g in 
the weak-coupling limit, where the proportionality factor 
can be found to be —2 from the Heine-Stieltjes connec¬ 
tion for di = 1/2. For all other levels j ^ i, Zj a remains 
finite and the dominant order is O(g 0 ). 

A dominant order of 0(1/g) in A, then corresponds 
to an excited state i in the noninteracting limit, while a 
dominant order of O(g 0 ) results in an non-excited state i. 
This behaviour can be generalized through a connection 
of A i to occupation numbers, as will be shown in the 
section about form factors. Note that the divergence in 
g —» 0 can pose numerical problems, so this suggests the 
use of gAi as variables instead of A,. 


C. Solving the equations 

Due to the similarity of Eq. ( f22| ) to the set of equa¬ 
tions found for the rational model”, it is straightforward 
to extend the solution method for the rational model to 
our equations. General sets of nonlinear equations have 
to be solved by an iterative approach starting from an 
initial guess, such as the Newton-Raphson method. This 
method converges quadratically to the solution if the ini¬ 
tial guess lies in the basin of attraction. An efficient 
numerical approach can be implemented based on this 
method once we have access to a sufficiently good initial 
guess for the solution. 

As in [I], we start from an approximation of the solu¬ 
tion in the weak-coupling limit (|g| <C). A solution at any 
value of the coupling constant can then be found by adia¬ 
batically varying g starting from the weak-coupling limit 
and using the solution at the previous step as the starting 
point for an iterative solution at the current step. This 
initial guess can be improved by using a Taylor expansion 
of the solutions at the previous step, 


(t \ • 

K(g + 6g)~ Ai{g) + 5g-^- (29) 

since the derivatives of the A^ can be found by solving a 
linear system once the set of A i is known, similar to the 
procedure followed in |T]. The equations for the deriva¬ 
tives can be found by deriving the set of equations to g 7 
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leading to 


A,: 


dK 

dg 


A, 


15A, 1 " 

S % 2 ^ J V % 


SA,; 

as 


(30) 


By taking the higher derivatives of the original set of 
equations, linear equations can be found for higher or¬ 
der derivatives of A , t , which allows a Taylor expansion up 
to arbitrary order. For an efficient numerical implemen¬ 
tation, combining the Newton-Raphson method with a 
Taylor approximation up to first order already offers a 
remarkable increase in speed. 


D. Inverting the transformation 


Although many properties of the XXZ systems can be 
written in terms of the eigenvalue-based variables Aj, as 
will be shown in the following section, it is still neces¬ 
sary to obtain the rapidities from these variables for the 
evaluation of e.g. certain form factors. For the rational 
model, Slavnov’s determinant formulsP^ is the building 
block of many such expressions, and this formula relies 
on the construction of a matrix expressed in the rapidi¬ 
ties. If we wish to obtain a general formalism, it should 
be possible to obtain the rapidities from the eigenvalue- 
based variables for any XXZ model. 

In order to be as general as possible, we propose a 
method to obtain the rapidities that does not rely on the 
explicit expression for X and Z except in the very last 
step. We introduce an auxiliary index r, corresponding 
to an auxiliary single-particle level or an additional un¬ 
coupled rapidity e r , such that the Gaudin algebra can be 
extended. Following Ortiz et alPand Eq. (211, each Z ia 
can be written as 


Z r 'j Z r 


(31) 


Once an explicit expression for Z tJ is known, Z ri = 
Z(e r ,ei ) can easily be calculated after choosing e r . In¬ 
stead of solving for the rapidities {x a } 1 it is now possi¬ 
ble to find a transformation that allows us to determine 
the set of Gaudin algebra elements {Z ra }. If these are 
known, the rapidities can always be found by solving each 
separate equation Z ra = Z(e r ,x a ) for x a . 

This expression can be used to show that 


N 


N 


A i — ^ ( Zi a — ^ ( 


r t z r iZ r 

7—7 

xJrr 


N 


= -NZ r . 


(T + Z 2 ri )Y 


Zri Z r 


(32) 


We can now define a polynomial with the full set (a = 
1... N) of Z ra as roots 


N N 

P{Z) = (Z - Z ra ) = Y PN-mZ m , (33) 

ol— 1 m—0 


with Pq = 1. Once the coefficients Pn-tu are known, 
the roots of this polynomial can be determined to find 
the variables. This method arises naturally for differ¬ 
ent problems in the theory of integrable systems, such as 
the Heine-Stieltjes connection^ the numerical meth¬ 
ods by Guan et alP 1 and Rombouts et alP^, and the 
weak-coupling limit in RG models 6 . 

The definition of P(z) can now be used to consider 


P'(z) 

P(z) 


J2m=0 m PN-m,Z m 1 _ _ 1 

Em=0 P N-mZ m Z ~ ' 


(34) 


Evaluating this equality in z = Z r i , we find that 

A» = - NZ ri + (r + Z 2 ri ) (35) 

resulting in 


P , (Z ri ) _ 1 _ A i + NZi 

P ( Z r i ) Zri — Z ra r + za 


The coefficients of the polynomial P(z) can then be found 
by solving a linear problem 


N-l 


Y PN-m [mPZ™- 1 - A t Z™ + (m- N)Z™ +1 ] 


m—0 


= A- NTZ*- 1 . 


(37) 


We obtain n equations for N < n variables, but if the 


set of A i variables can be written as Eq. (32), these are 


linearly dependent and we are free to choose N equations 
from the full set and solve these. Once these coefficients 
are known, the variables Z ra can be determined using a 
root-finding algorithm such as Laguerre’s method. Al¬ 
though easy to implement, this method has the disad¬ 
vantage of being prone to numerical errors. Indeed, it 
is well-known that the roots of a polynomial are highly 
sensitive to changes in the coefficients, sometimes even 
changing the solutions at the qualitative level (a set of 
complex conjugate roots can be found numerically in¬ 
stead of two separate real roots). These considerations 
are not pressing for a limited number of excitations, but 
become more and more important for an increasing num¬ 
ber of excitations. 

This problem was also encountered for the ratio¬ 
nal rnodel^, after which an improved method was pro¬ 
posed by decomposing P(z) in the basis of Lagrange 
polynomial^!. This led to an increased accuracy and al¬ 
lowed the method to tackle problems with a large num¬ 
ber of excitations (up to a few hundreds). It is expected 
that these results can be generalized towards our prob¬ 
lem because of the similarity of all necessary equations. 
We refer the reader to [2j for a detailed analysis of the 
method. 

As an illustration, the eigenvalue-based variables and 
the related rapidities have been plotted for the different 
models in Figure[l]as a function of the coupling constant. 
For each model, the levels c, are given by a picket-fence 
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modef^ with equal level spacing: ej = For each 

system singular points occur, where multiple real rapidi¬ 
ties combine and continue as a pair of complex conjugate 
variables. As is clear from Figure [l] the RG equations 
(Eq. [l2] ) become highly singular, whereas the eigenvalue- 
based variables vary very smoothly. 


E. Degenerate models 


The discussion has been limited to di = 1/2 models 
only thus far. This is sufficient for the majority of inter¬ 
acting (quasi-)spin systems, however in situations with 
higher symmetries arbitrary degeneracies di > 1/2 may 
occur. For the rational model, the set of equations for 
the eigenvalue-based variables were derived starting from 
a differential equation, where each equation for A, is ob¬ 
tained by evaluating this equation at a different level e t . 
For degenerate levels, additional equations could be ob¬ 
tained by taking derivatives of the differential equation 
and evaluating these at the different values ep-A This 
necessitated the introduction of new variables 



N 


Ai X a ) n 


(38) 


which are highly reminiscent of the set of variables in¬ 
troduced by Rombouts et alP^. The total number of 
variables is determined by the number of single-particle 
levels weighted by their degeneracies. For each single 
particle level i with degeneracy 2 di + 1, taking the first 
2 di — 1 derivatives of the differential equation and eval¬ 
uating each equation at ej then results in a closed set of 
equations for the set of variables Aj, A- 2 ^,..., A- 2di \ V*. 

The outlined procedure is slightly more involved for 
the class of XXZ RG models, however it is analogous to 
the rational case, and remains therefore tractable. The 
main idea is identical: starting from a continuous rep¬ 
resentation of the equations it is possible to obtain any 
number of equations for the total set of variables. The 
continuous representation of the variables is given by 


for A i if di = 1 / 2 , but for larger degeneracies the set of 
equations also depends on 





a=l 


(41) 


By taking the derivative of Eq. (401, additional equations 
can be obtained linking these higher-order variables to 
the original variables. The derivation of this method is 
given in Appendix [B] It is remarkable that the derivative 
of Mz), a summation over Z(z,x a ), can still be related 
to A^ 2 )(z), a summation over Z(z,x a ) 2 . 

As an example, the evaluation of the first derivative 
results in 


Aj (AT + A| 2) ) = - X - (AT + Af 

+ E - A?) 

+ E djZij (Af } + r) 

+r(l — di)Ai + (1 — Gp)A- 3) . (42) 


For di = 1, we have obtained a closed set of equations 
in {A^Ap*} for each level. For arbitrary di, the first 
2 di — 1 derivatives results in a closed set of equations. An 
additional equation for the total number of excitations 
can easily be determined as 


N = ~Y J 9diK (43) 

i —1 


As an illustration, the results for a level-independent re¬ 
duced BCS model describing neutron superfluidity 18 in 
56 Fe are given in Figure [ 2 J The energy levels are deter¬ 
mined as the levels of a spherical symmetrical Woods- 
Saxon potential and the di, dependent on the angular 
momentum quantum numbers, vary from 1/2 to 5/2. 


IV. FORM FACTORS 


N 

H z ) = E Z ( z ’ x <*)’ ( 39 ) 

a=l 


where A(e^) = Aj. A similar derivation as for the di = 
1 /2 model in Appendix[A]results in a continuous equation 


A (~)] 2 = TN ^1 - N + 2 ^ djj - -A (z) 

+ 2^2 d j z (z,e j ){A( z ) -A (ej)) 

+ ^2 Z(z, x a ) [Z(z, x a ) - 2diZ(a, x a )], (40) 


which holds if z / ej for all j 7 ^ i. The evaluation of these 
equations at z = £i results in the known set of equations 


A. Overlap with Slater determinants 

When performing calculations with Bethe Ansatz 
states, it is customary to expand them in a basis set 
of uncorrelated wavefunctions. These expansion coeffi¬ 
cients are given by permanents of matrices in the general 
case, following from the theory of multivariate generating 
functiontPHl. However, these expressions are not practical 
for computational purposes since the evaluation of a per¬ 
manent scales exponentially with the matrix size. Com¬ 
pared to the determinant, which can be evaluated in a 
polynomial time, permanents require a factorial scaling 
computational time and therefore severely limit the size 
of the systems that can be considered. 

Luckily, multiple formulae exist for the rational model 
linking the permanents to determinants, allowing numeri¬ 
cally efficient expansions of the Bethe Ansatz state. Here, 
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Figure 1. Evolution of the variables for the ground state of a ’picket-fence’ model with 6 excitations in 12 doubly-degenerate 
levels [1., 2., 3.,..., 12.]. The rational, the hyperbolic and the trigonometric models are organized in columns, while the rows 
depict the set of eigenvalue-based variables gAi and the real and imaginary part of the rapidities as a function of the coupling 
constant g = 0 • • • — 1. The used parametrizations are given by Eqs. H>, (§ and Q respectively. Note the Moore-Read point 
at g = — 2/n = —0.167 for the hyperbolic modeP, where all rapidities condense to zero, and the eigenvalue-based variables all 
become equal to —1. Due to the periodicity of the trigonometric model, the real parts of the rapidities lie within the interval 
[— 7 t / 2 , + 7 t / 2 ], as marked in the figure. 


we present several such results for di = 1/2 XXZ systems 
as a generalization of the results previously obtained for 
the XXX modeP and the inhomogeneous Dicke modeP. 
Starting from these expressions, determinant expressions 
for systems with arbitrary di can also be obtained as a 
limiting case where several rapidities coincide, but this 
will not be discussed here. 

Multiple determinant formulae exist for the overlap of 
a Bethe Ansatz state in the rational model with a Slater 
determinant 


N 

no•••+>> ( 44 ) 

a—1 


which is an uncorrelated wavefunction and an eigenstate 
in the uncoupled limit (g = 0). The Bethe Ansatz state 


([TT]) can be expanded in this basis as 


N / n \ 

i^> = n i 0 ) 

OL—l \ i—1 ) 


= £|{U) <{UhM 

{ia} 

(45) 

^ . |{G}) per (Xy aCK ) , 

{U 

(46) 

where {i a } runs over the partitioning of n levels 
different excitations modes and 

over N 

({*a}|?M = per i\X ia0 \), 

(47) 


which is the permanent of an (N x N) matrix with matrix 
elements given by Xi a0i . This is a formula independent 
of the integrability or the explicit expression for X and 
is simply a result of the structure of the wavefunction. 
However, by imposing the Gaudin algebra on the matrix 
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Figure 2. Evolution of the variables for the ground state con¬ 
taining 11 excitations in 10 levels for a varying (negative) cou¬ 
pling constant g, describing the pairing of neutrons in 56 Fe by 
means of a level-independent reduced BCS model. The pa¬ 
rameters of the system are taken from [18], where the set 
of di are given by [3/2,1,1/2, 2,1, 3/2,1/2,1/2,1, 5/2] with a 
maximal degeneracy of 2 di + 1 = 6. The first figure shows 
the evolution of the eigenvalue-based variables, while the sec¬ 
ond and third figure correspond to the real respectively the 
imaginary part of the rapidities. 


elements, it is always possible to rewrite this permanent 
as a determinant. 


per([X ioa ]) 


(ria> 6 ^+*i>) 

ria.a Xi a a 


det [ X iaa * X ia „](48) 


where the Hadamard product of two matrices is intro¬ 
duced, defined as (A * B) b j = AijBij. This is a gen¬ 
eralization of the Borchardt or Izergin determinant rep¬ 
resentation for the rational rnodePB®] an d has a simi¬ 
lar structure to some results presented previously for the 
trigonometric rnodeP^. 

In ref. [3], a case was made for a RG theory that would 
require only the eigenvalue-based variables instead of the 
rapidities. This would be desirable from a numerical 


point of view as the eigenvalue-based variables are free of 
singularities, opposed to the singularity-prone rapidities. 
In addition, this would enable us to skip the inversion 
step, which is the main bottleneck of this method. 

For the rational model, it was shown by Faribault and 
Schuricht® that 


per 



= det J 


(49) 


where the left-hand size is the permanent of an TV x N 
Cauchy matrix and the right hand side is the determinant 
of a matrix J defined as 


Jab — 


t-^OL = 


a=l Ci 
1 


c-E 

. —x n ' 


N 


c^a 6i a —ei i 


if a = b 
if a/i) 


(50) 


where the only dependence on the rapidities is through 
the diagonal elements, which can be expressed in terms 
of the eigenvalue-based variables. The ingenious proof 
of this expression is based on a recursion relation which 
was found to hold for both sides of Eq. (501. We will 


show that these results can be easily generalized to XXZ 
models. 

The Gaudin equations can be used to rewrite the per¬ 
manent in the XXZ model in a structure similar to a 
Cauchy matrix. We once again introduce an auxiliary 
level r and use the Gaudin equations to write 




A-ri A- rex 

Z r i Z rr . 


r/ r + Z r i Z ra , . 

Z ia = y _ 7 -■ (51) 

7 rr * 


Rewriting all matrix elements in the permanent <fi allows 
us to take 


per ([W a a]) = per 


Z-ria X T c* 
7 7 



resulting in a Cauchy matrix, which we can rewrite as a 
determinant 


per {[Xi a a]) 



(53) 


with J redefined as 

T _ j E«=i z—- Ez—h—■ *« = & 
ab \ z-h- if 

( 54 ) 

By multiplying each row and column c with X r i c and 
compensating for these factors in the prefactor, this can 
be written as 


per ([Xi aa ]) 


EL Xroc 

ria E+ 


det J 


(55) 


with 


Jab — 


^iV xl ia 

Z ri -Z rci 

X r i x riu 


E 


N Xl ia 
C^a Z r i a —Z r i c 


if a = b 
if o^6, 
(56) 
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where we recognize X v j (j in the off-diagonal elements 
and the diagonal elements can be rewritten by using the 
Gaudin algebra until 

T f E a =l ^iaa ~ Kc + Z r i a if a = b 
ab \x ia i b if a±b' 

(57) 

The only explicit dependency on the rapidities is now 
found in the prefactor and in the diagonal elements, 
which only depend on the eigenvalue-based variables. 
Fortunately, the prefactor can be absorbed in the defi¬ 
nition of the state without loss of generality, which will 
still allow the calculation of all necessary normalizations 
and prefactors. Note the appearance of in the diag¬ 
onal elements, which can be easily evaluated and links 
the matrix to the prefactor for the Bethe Ansatz state 
through our choice of the additional level r. 

To recapitulate, the overlap of a state 

!{*«»= nfzfM !+...+> (58) 

a =1 \i=l ^ ra ) 

with a Slater determinant 


The RG equations for the dual state are given by 


n-N 

- 1 + | Z Ja f = g ^2 z 0'a ', Va' = 1 ...n-N 

i 

(64) 

or, written in the eigenvalue-based variables, 

[A'] 2 = iV(n—iV)r+-A'+^ Zji(. A'- A'), Vi = 1... n. 

" 3+i 

( 65 ) 

Both representations describe, up to the normalization, 
the same state. So by comparing the eigenvalues of the 
constants of motion 


n = 


1 

2 


1 — gAi + 



1 

2 


1 - f/A' + 7:^2 Zik 

Z k^i 


( 66 ) 


it can be seen that the dual representation of an eigen¬ 
state is determined by a set of dual eigenvalue-based vari¬ 
ables given by 


KU> 


N 


k*i ■■•**}>=m 

3= 1 


is given by 


({ii.. .ijv}|K}) 


det J 
Ha Xri 


with 


A ia - Z iaic + Z ria if a = b 
X% a i b if cl ^ b 


(59) 



Vi = 1... n, 


(67) 


which can be verified by simply substituting these vari¬ 
ables in the dual eigenvalue-based equations. The exis¬ 
tence of the dual state allows one to write the overlap 
of a normal state and a dual state as the overlap of a 
Bethe Ansatz state with the fully-filled state. This will 
be exploited in the following section to obtain the nor¬ 
malization of the Bethe Ansatz states and several form 
factors. 


(61) 

C. Normalization and form factors 


B. Dual representations 


So far, all Bethe Ansatz states have been created by 
acting with creation operators on the vacuum state, de¬ 
stroyed by all annihilation operators. Due to the particle- 
hole symmetry of the RG models, every Bethe Ansatz 
state can also be constructed by acting on the fully-filled 
state, annihilated by all creation operators, with a set 
of generalized annihilation operators, the so-called dual 
representation. 

A renormalized Bethe Ansatz state is given in its nor¬ 
mal and dual representation as 


IK » 

IK'}) 


nlEfM IKK, ( 62 ) 

a=l \i= 1 ra J 

n-N / n Y \ 

n ) it - - -1>. ( 63 ) 

a' = 1 \i=l ra J 


Due to the dual representation, we have two different 
representations of each Bethe Ansatz state. These two 
states have different normalizations and can be written 
as 

Sjl \i---i)=N a \{x a }) n 


n-N / 

ik» = n 

a' — 1 \ 

where \{x a }) n is the normalized eigenstate. This can be 
used to calculate 


X, ; , 


X r 

i=l r 


-Si It t ) = x a . | K» r 


N 


ik»= n e 


CK—1 \ i—1 


x i: 

X r , 


(K'llKI) = N a N a > 

n-N / n v \ N / n v \ 


— l \ i—1 


a—1 \i=l 


det J 

TEk 


( 68 ) 
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with 


2A i + | — 'Yhk^i Z ik + Zri if * — J 
Xij 


(69) 


since the overlap between the two states is just the 
overlap of a Bethe Ansatz state determined by {x} = 
{x Q }U{x a -} or by A * = Aj+A' = 2Ai+2/g with |f ■ • • t)- 
The ratio N a /N a i can easily be found by taking the over¬ 
lap of both representations with a single reference state 
and taking this ratio. Knowing both the ratio and the 
product of the normalizations, the normalization of both 
representations is uniquely defined 

In a similar manner as for the XXX modeP®, these 
results can be used to calculate form factors. Since this is 
a generalization of the results presented in these papers, 
we will only summarize the results for the XXZ model 
here. 

The form factor for the raising operator between a 
state {x Q } with N excitations and a dual state {x(,} with 
N + 1 excitations is given by 

, det 

= n y (70) 

II i^k ^ri 

with 


jk _ 
J ab ~ 


_ J A“ + A f + g ~ YZx,ijtk Zil + z ri if O — 6 
X a b if a / b 

( 71 ) 

and a, 6 ^ k. The form factor for a lowering operator is 
then given by the Hermitian conjugate of this expression. 

The form factor for local operators 5° can similarly be 
found, both for diagonal and off-diagonal expectation val¬ 
ues. Here we will make the dependency on the coupling 
constant explicit and define {£ a (g)} as an eigenstate at 
coupling constant g. 

From the Hellmann-Feynman theorem we obtain for 
the diagonal expectation values 


(W a (9)}\s° k \{x a m = 

1 / f)\ a \ 

2{~ 1 + g 2 ^i) 


(72) 


while the off-diagonal expectation values can be found as 

({<mS° k \{x a (g)}) 


= l -{gAl-gKl + 2) 


(WJg)}\{x a (g + dg)}) 


dg 


1 _ n _ f) 

= -( 5 A^- 5 A^+2)^^AdetJ fc 

k= 1 


dg 


(73) 


with 



-E 


l^i 


if a = b 
if a b 


(74) 


and a,b ^ k. These form factors can be written in an 
expression only dependent on the eigenvalue-based vari¬ 
ables, so the rapidities do not need to be determined 


explicitly. However, if these have been determined, other 
form factors can also be calculated by commuting the 
operator through the creation operators in the Bethe 
Ansatz state, as has been shown in [35]. This allows 
the computation of form factors such as e.g. S' k Si as a 
sum over determinants. 

It is interesting to note that the expectation value for 
Sf offers a physical interpretation for Aj, since the deriva¬ 
tive of Aj to g fully determines the filling of the level i 
as described by {Sf). Knowledge of the evolution of Aj 
with a changing coupling constant is then equivalent to 
knowing how the N excitations are distributed over the 
n levels. 


V. CONCLUSIONS AND OUTLOOK 

In the present paper, we presented a numerical solution 
method for the RG equations for XXZ integrable models, 
while also proposing determinant expressions for the nor¬ 
malization and form factors of the Bethe Ansatz states 
that are only dependent on the new set of eigenvalue- 
based variables. The availability of an efficient solution 
method and explicit expressions for the overlaps opens up 
possibilities for the investigation of integrable and non- 
integrable quantum systems. 

An efficient numerical expression for overlaps be¬ 
tween two Bethe Ansatz states allows for a numeri¬ 
cally nearly exact investigation of the dynamics result¬ 
ing from a quantum quench in the reduced BCS pairing 
modepS while for models close to integrability the eigen¬ 
states of the integrable model are being used to study 
decoherenc^H. These results were all obtained for the 
rational XXX model, and it would be interesting to in¬ 
vestigate similar dynamics for XXZ models, where a rich 
phase diagram is knowrP. 

Due to the favourable computational scaling of this 
method it is also possible to numerically approach the 
thermodynamic limit and investigate the limiting be¬ 
haviour of the system. Such an investigation has been 
carried out for the reduced BCS Hamiltonian by El Araby 
and BaeriswyP^, where exact results for large system 
sizes were compared with the BCS Ansatz. Again, it 
would be worthwhile to compare the results rigourously 
with the results from BCS mean-field theory®. 

The overlaps between Bethe Ansatz states and Slater 
determinants have also been used in quantum chemistry 
for strongly correlated systems. These can not be accu¬ 
rately described by means of a single Slater determinant, 
and several wavefunctions have been proposed to cap¬ 
ture the correlation present in the system. One example 
is the class of geminal-based wavefunctions, which have 
a structure highly reminiscent of the Bethe Ansatz states 
found in the Gaudin models. Unfortunately, these wave- 
functions are only computationally feasible if a numer¬ 
ical efficient expression is known for the overlaps with 
Slater determinants. One such class (APr2G) is based 
on the rational XXX modeP^SJ, and the presented re¬ 
sults should allow for an extension of these geminal-based 
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wavefunctions with a class based on the XXZ models. 
The results presented here should also allow for a varia¬ 
tional approach based on the Bethe Ansatz states in the 
XXZ model. By starting from Richardson’s general so¬ 
lution for the Gaudin algebra (Eq. [8]) and varying over 
the free parameters in the wavefunction, it is possible to 
scan over all possible Gaudin models (XXX and XXZ), 
allowing more freedom and as such a better approxima¬ 
tion to the energy compared to a variational approach 
based solely on the XXX model. 

It was already noted by Babelon and Talalaev that 
similar equations could be found for the XXX Heisen¬ 
berg spin chairP^l, so the question naturally arises if 
this method can be generalized to other Bethe Ansatz 
equations obtained from the quantum inverse scattering 
method^. 


which can again be rewritten by making use of Eq. 21 


A? =-N(N - 1)T + - y^Z ai 
9 a 

+ EE (r + Zji(Zja + z ai )) 

a j^i 

= N(n - N)T - -Aj + 'Y ZjiiAj - A;), (A3) 
9 fri 


resulting in a set of equations closed in the variables 
A i, i = 1... n. 


Appendix B: Obtaining equations for degenerate 
models 
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Here, we show how to obtain additional equations for 
degenerate models starting from a continuous equation 
for the eigenvalue-based variables. We define 

N 

A, (P) = Y Z ^ (Bl) 

a=l 


Appendix A: Deriving the equations 

In this appendix we derive the equations for 
the eigenvalue-based variables by making use of the 
Richardson-Gaudin equations and the Gaudin algebra. 
All summations with indices labelled by Greek charac¬ 
ters run over the N rapidities, while for indices labeled 
by Latin characters the summations run over the n en¬ 
ergy levels. Starting from the definition of the variable 
Aj, we can write 

A i = 'y ] Zi a Zi[) = y ] Zi a Zip + ^ ' Z ia 

ol, fi ol,(3^ol ol 

= — Yj [r + z a p(z ai + z ip )] + Y z l a , (Ai) 

ol,(3^ol ol 


where we have used Eq. 21 to rewrite Z{ a Zip. The anti¬ 
symmetry of Zi^ and Z a p can be used to obtain 


A 2 i=- Y (r + 2 z a pZ ai ) + Yz, 

a.,f3^0i ol 


= —N(N -l)T + 2Yz m [Y Z ^} +Y Z 


Making use of the RG equations (Eq. 12), we obtain 


Ai = —N(N — 1)T + 2 22 Zai ( ~ g + \Y Z i« | + £3 


= —N(N — 1)T H— Y Zai +YJ2 ZaiZ i a ’ 

Q ol OL j^i 


and 

N 

A to(z) = Y z ( z > x c') p , pe N. (B2) 

a —1 


By taking the derivative with respect to z of Eq. (401, 
additional equations can be found. Unfortunately, an 
explicit expression for Z(z,6i) is needed if we want to 
relate the derivative of A (z) to A^ 2 ) (z) and obtain a closed 
set of equations. In order to circumvent this problem, we 
introduce a fixed auxiliary level e r , so that every Z can 
now again be parametrized as in Eq. (|2l|), 


Z(z,x a ) 


Z(e r ,z)Z(e r ,x a ) +T 
Z(e r ,z) - Z(e r ,x a ) 


(B3) 


so instead of deriving this equation with respect to the 
continuous variable z, we can derive with respect to 
Z(e r ,z) and multiply with — A 2 Z , resulting in an equa¬ 
tion independent of the explicit expression for Z. For 
notational ease, we take Z rz = Z(e r ,z). Then 


-X 2 z -^—Z(z,x a ) = X(z,x a ) 2 = Z(z,x a ) 2 + T, (B4) 

dZj ^2 

which is independent of the auxiliary level. The action 
of this operator (derivation and multiplication) is given 

by 


A (p >(*) -a P A (p+1 \z) + pTA^-^iz), (B5) 


and 


(A2) 


Z{z,e 3 ) p ^pZ{z,e^ p+1) + pTZ(z, e i ) (p - 1) . (B6) 
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By acting with this operator on the equations for A(z), 
we obtain 


A(z)(AT+A^(z)) = -- (NT + A (2 >(z)) 


+ + Z( z > e j) 2 ) (A(z) - A(e,-)) 

+ y( a ( 2 ) (z) + r) 

T ^ ^ *£a) i^( z i ^a) djZ (tzi . X a )) 

a 

+r(A(z)-d i A(e i )) (B7) 


resulting in Eq. (42) when evaluated in z = e^. Succes¬ 
sive derivatives result in as many additional equations as 
there are variables, allowing a closed set of equations to 
be obtained. 
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